#include "./crystal.h"
#include <math.h>


static
Real mds_crystal_meshes( Real x )
{
	x = mds_anint(x);
	if( x <= 1) x = 1;
	return x;
}

size_t mds_crystal_close_packed_count( mds_uvector *cr, const size_t n, const mds_vector *box )
{
	assert( n > 0 );
	assert( box != NULL);
	assert( box->x > 0 );
	assert( box->y > 0 );
	assert( box->z > 0 );
	assert( cr != NULL );
	{
		const Real d     = (box->x*box->y*box->z) / ( n );
		const Real alpha = (Real)pow( MDS_F(1.0)/(MDS_F(4.0)*d), MDS_F(1.0)/MDS_F(3.0) );
		
		/*
		 fprintf( stderr, "d = %g\n", d);
		 fprintf( stderr, "alpha= %g for [%g %g %g]\n", alpha, box->x, box->y, box->z);
		 */
		
		cr->x = mds_crystal_meshes( alpha * box->x);
		cr->y = mds_crystal_meshes( alpha * box->y);
		cr->z = mds_crystal_meshes( alpha * box->z);
		
		return cr->x * cr->y * cr->z * 4;
	}
}



void mds_crystal_build_fcc( mds_vector_set r, const size_t n, const mds_uvector *cr, const mds_vector *box, Real noise_level )
{
	assert( r  != NULL );
	assert( n > 0 );
	assert( cr != NULL );
	assert( box != NULL );
	assert( cr->x * cr->y * cr->z * 4 == n );
	{
		const Real dx = box->x / cr->x; const Real dx2 = MDS_F(0.5) * dx;
		const Real dy = box->y / cr->y; const Real dy2 = MDS_F(0.5) * dy;
		const Real dz = box->z / cr->z; const Real dz2 = MDS_F(0.5) * dz;
		const mds_vector noise = { noise_level * dx, noise_level * dy , noise_level * dz, {0} };
		mds_vector u;
		size_t     m = 0;
		size_t     i,j,k,iref;
		/** SUBLATTICE A */
		
		r[0].x = 0;
		r[0].y = 0;
		r[0].z = 0;
		
		/** SUBLATTICE B **/
		
        r[1].x =  dx2;
        r[1].y =  dy2;
        r[1].z =  0;
		
		
		/** SUBLATTICE C **/
		
        r[2].x =  0;
        r[2].y =  dy2;
        r[2].z =  dz2;
		
		
		/** SUBLATTICE D **/
		
        r[3].x =  dx2;
		r[3].y =  0;
        r[3].z =  dz2;
		
		
		/** Repeat sublattices */
		for( k=0; k < cr->z; ++k )
		{
			const Real z = k * dz;
			for( j=0; j < cr->y; ++j )
			{
				const Real y = j * dy;
				for(  i=0; i < cr->x;  ++i )
				{
					const Real x = i * dx;
					for(  iref=0; iref < 4; ++ iref )
					{
						const size_t p = iref+m;
						assert( p < n );
						r[p].x = r[iref].x + x;
						r[p].y = r[iref].y + y;
						r[p].z = r[iref].z + z;
					}
					m += 4;
				}
			}
		}
		assert( m == n );
		
		/** create noise */
		for( i=0; i < n; ++i )
		{
			mds_vector *r_i = r+i;
			mds_vector_unit( &u );
			r_i->x += u.x * noise.x;
			r_i->y += u.y * noise.y;
			r_i->z += u.z * noise.z;
			
			mds_pbc_xyz( r_i, box);
		}
	}
	
}

